function lower_bound = gmmBayesLowerboundLosses(mix, data, vars)


K = mix.ncentres;

% E[ln p(X|Z, mu, Lamba)]
q1 = zeros(1, K);
for k = 1:K
    vars2 = vars;
    vars2.Z(:, k) = 0;
    vars2.Z = vars2.Z./(sum(vars2.Z, 2)*ones(1, K));

    [mix.sample_means mix.sample_covars] = computeSufficientStatistics(data, vars2, mix.covar_type, 1e-6);
    
    q1(k) = ElnpX_ZmuLambda(mix, vars2);
end

% E[ln p(Z|pi)]
q2 = ElnpZ_pi_losses(mix, vars);

% E[ln p(pi)]
q3 = Elnppi_losses(mix);

% E[ln p(mu, Lambda)]
q4 = ElnpmuLambda_losses(mix);

% E[ln q(Z)]
q5 = ElnqZ_losses(vars);

% E[ln q(pi)]
q6 = Elnqpi_losses(mix);

% E[ln q(mu, Lambda)]
q7 = ElnqmuLambda_losses(mix);


lower_bound = q1 + q2 + q3 + q4 - q5 - q6 - q7;



% lower_bound
% 
% [mix2 vars2] = killComponents(mix, vars, 5);
% ll = gmmBayesLowerbound(mix2, data, vars2, 1e-6)
% 
% [mix2 vars2] = killComponents(mix, vars, 3);
% ll = gmmBayesLowerbound(mix2, data, vars2, 1e-6)
% 
% pause
